Visualización y geolocalización de datos

Máster en Data Science & Business Analytics (con R Software)

Objetivos de aprendizaje

¿Por dónde empezamos? Recursos interesantes

Libros de referencia:

• Spatial Data Science with applications in R

• Geocomputation with R

• Displaying time series, spatial and space-time data with R

Recursos de estadísitca espacial en R:

• rspatial
• R-spatial

Otros recursos web interesantes:

• Libro de Ruben
• Rositter
• Libro Hengl
• Libro de Páez

1 La revolución de los geodatos

Que estamos en la la era del dato, que los datos son el petroleo del siglo XXI y que estamos rodeados de datos es una cuestión que ya hemos hecho inherente a nosotros. Estamos en el momento del dato, donde la profesión de Data Scientist se ha convertido en la profesión más sexy del siglo XXI según vaticinó en 2012 Harvard Business Review. Cada segundo se producen 1,7 MB de datos/persona y cada año esta cifra se duplica se duplica.

Este incremento exponencial de los datos ha sido posible, sin duda, gracias al desarrollo de la tecnología, los ordenadores, los teléfonos móviles móviles, los satélites, internet, etc… y asociado a estas nuevas herramientas, se ha producido una lluvia sin precedentes hasta el momento de datos espaciales o datos georreferenciados. Cada teléfono inteligente tiene un receptor de posicionamiento global (GPS) y una multitud de sensores en dispositivos que van desde satélites y vehículos semiautónomos hasta científicos ciudadanos que miden incesantemente cada parte del mundo. La tasa de datos producidos es abrumadora. Un vehículo autónomo, por ejemplo, puede generar 100 GB de datos por día (The Economist, 2016). Los datos de teledetección de satélites se han vuelto demasiado grandes para analizar los datos correspondientes con una sola computadora.

Esta revolución de los geodatos y el análisis de los datos espaciales no sólo impulsa la demanda de hardware informático de alto rendimiento y software escalable y eficiente para manejar y extraer la señal del ruido, lo que se conoce como Geocomputación, sino que ha dado lugar una nueva rama de conocimiento, Spatial Data Scicene (SDS) o Ciencia de Datos Espaciales.

Como ejemplo, el Ministerio de Transportes, Movilidad y Agenda Urbana llevó a cabo durante los años 2020 y 2021 el denominado Estudio de movilidad con Big Data, cuya fuente principal de datos fue el posicionamiento de los teléfonos móviles anonimizado. Estos datos permiten por ejemplo analizar la movilidad entr diversas zonas del territorio español de manera diaria:

2 Geocomputación y R-spatial

2.1 Geocomputación. ¿Qué es? PUEDO EMPEZAR YO Y LUEGO DIEGO LE DAS UNA VUELTA. ver lo que tiene DIEGO

La geocomputación es un término relativamente nuevo pero influenciado por otros términos clasiscos. Es decir, puede ser visto como una parte de la Geografía, la cual tiene más de 2000 años de historia (Talbert_2014) o como una extensión de los Sistemas de Información Geográfica comunmente conocidos por GIS (del inglés, Geographic Information Systems) (Neteler and Mitasova 2008), los cuales emergieron en la década de los 60 (Coppock and Rhind 1991).

Según Lovelace, Nowosad, & Muenchow (2019), el término Geocomputación se utilizó por primera vez en 1996 en la primera conferencia que trataba de Geocomputación. Abrahart, Openshaw, Abrahart, & See (2000) definenen la geocomputación como una herramienta que trata de utilizar los diferentes tipos de datos geográficos y de desarrollar herramientas geográficas relevantes dentro del contexto general de un enfoque ‘científico.’

La geocomputación está estrechamente relacionada con otros términos que incluyen:

• la ciencia de la información geográfica (GIScience),

• la geomática,

• la Geoinformática,

• la ciencia de la información espacial

• la ingeniería de Geoinformación (P. Longley 2015)

• y ciencia de datos geográficos (GDS, Geographic Data Science).

Cada término comparte un énfasis en un enfoque científico (que implica reproducible y falsable) influenciado por los SIG, aunque sus orígenes y principales campos de aplicación difieren.

La ciencia de datos geográficos, por ejemplo, enfatiza las habilidades de “ciencia de datos” y grandes conjuntos de datos, mientras que la geoinformática tiende a enfocarse en estructuras de datos. Pero las superposiciones entre los términos son más grandes que las diferencias entre ellos y usamos geocomputación como un sinónimo aproximado que los encapsula a todos: todos buscan usar datos geográficos para trabajos científicos aplicados.

Resumir del libro Lovelace et al. (2019). Empezar ilustrando con un ejemplo, O HACER REFERENICA AL DE LOS AVIONES

2.2 El ecosistema R-spatial y datos geográficos en R o algo así…

3 Estadística espacial

(lo tengo literal del artículo mio de AHEPE, resumir)

La estadística espacial reconoce y aprovecha la ubicación espacial de los datos a la hora de diseñar, recopilar, gestionar, analizar y mostrar las observaciones. Éstas son generalmente dependientes, si bien existen modelos espaciales a disposición del investigador que permiten tratar con dicha dependencia espacial a la hora de llevar a cabo labores de predicción. La estadística espacio-temporal incorpora, además, el tiempo y su interacción con el espacio como argumento de ayuda en tales labores predictivas.

Las mediciones y modelos espaciales están presentes, sorprendentemente, en una amplia variedad de disciplinas científicas. Los orígenes de la vida humana vinculan los estudios de la evolución de las galaxias, la estructura de las células biológicas y los patrones de asentamiento arqueológicos. Los ecologistas estudian las interacciones entre plantas y animales. Silvicultores y agricultores necesitan investigar las variaciones que se producen en el terreno para sus experimentos. La estimación de las precipitaciones y de las reservas de oro y petróleo es de vital importancia económica. Estos son, entre otros, buenos ejemplos de la importancia del espacio (espacio-tiempo en su caso) en el mundo de la Ciencia.

En todo caso, la geología, la edafología, el tratamiento de imágenes, la epidemiología, la agronomía, la ecología, la silvicultura, la astronomía, el estudio de la atmósfera, la economía, o simplemente, cualquier disciplina que trabaje con datos espaciales recopilados de diferentes lugares y en distintos instantes temporales, necesita del desarrollo de modelos geoestadísticos que indiquen la estructura e intensidad de la dependencia espacial y/o espacio-temporal presente en los fenómenos que comprenden.

Sin embargo, el estudio de la variabilidad espacial, y sobre todo espacio-temporal, es una disciplina relativamente nueva en el marco de la Estadística, lo que explica la escasez de instrumentos de estadística espacial 30 años atrás. En los últimos 10 años ha habido una creciente toma de conciencia de esta necesidad, habiéndose realizado un gran esfuerzo por buscar herramientas adecuadas y útiles a tales efectos. Y todo ello porque utilizar modelos espaciales o espacio-temporales para caracterizar y explotar la dependencia espacial (o espacio-temporal) de un conjunto de observaciones tiene importantes ventajas:

1. Modelos más generales, ya que, en la mayoría de los casos, los modelos clásicos que no tienen en consideración la dimensión espacial o la interacción de las dimensiones espacial y temporal son un caso particular de un modelo espacial o espacio-temporal.

2. Estimaciones más eficientes: de la tendencia, de los efectos de las variables explicativas, de promedios regionales,…

3. Mejora de las predicciones: más eficientes, con propiedades de extrapolación más estables,…

4. La variación espacial no explicada en la estructura de la media debe ser absorbida por la estructura del error, por lo que un modelo que incorpore la dependencia espacial puede decirse que está protegido frente a una mala especificación de este tipo. Esto, en muchos casos, tiene como resultado una simplificación en la especificación de la tendencia; en general, los modelos con dependencia espacial suelen tener una descripción más parsimoniosa (en ocasiones con muchos menos parámetros) que los clásicos modelos de superficie de tendencia.

Estas mejoras de la estadística espacial y espacio-temporal, junto con el fuerte y reciente desarrollo de los Sistemas de Información Geográfica o GIS (Geographic Information System), han propiciado que en la actualidad exista una importante motivación por la búsqueda de herramientas espaciales o espacio-temporales.

3.1 Antes de continuar… dependencia espacial.

Frecuentemente los datos tienen una componente espacial y/o temporal asociada a ellos y es de esperar que datos cercanos en el espacio o en el tiempo sean más semejantes que aquellos que están más alejados; en cuyo caso no deben ser modelados como estadísticamente independiente, sino que habrá que tomar en cuenta esa dependencia espacial o espacio-temporal.

De forma natural y de acuerdo a la Ley Tobler (1973) surge la idea de que los datos cercanos en el espacio o en el tiempo serán más similares y estarán más correlacionados entres sí que aquellos que están más lejanos. Además, esta correlación disminuye al aumentar la separación entre ellos, por lo que se puede pensar en la presencia de una dependencia espacial o espacio-temporal. Esto da lugar al concepto de proceso espacial o espacio-temporal.

Si los datos no exhiben dependencia espacial no tiene sentido aplicar las herramientas de estadística espacial. Veamos un ejemplo simulado de unos datos que muestras dependencia espacial y otros puramente aleatrorios.

library(geoR)
#> 'RandomFields' will use OMP
library(fields)

par(mfrow = c(1, 2))

set.seed(2022)
sim1 <- grf(441, cov.pars = c(1, 0.6))
#> grf: simulation(s) on randomly chosen locations with  441  points
#> grf: process with  1  covariance structure(s)
#> grf: nugget effect is: tausq= 0 
#> grf: covariance model 1 is: exponential(sigmasq=1, phi=0.6)
#> grf: decomposition algorithm used is:  cholesky 
#> grf: End of simulation procedure. Number of realizations: 1
points.geodata(sim1, main = "Dependencia espacial (positiva) (i)", col = tim.colors(), cex.max = 3)

# Independencia
set.seed(2022)
sim3 <- grf(201, cov.pars = c(0.01, 0))
#> grf: simulation(s) on randomly chosen locations with  201  points
#> grf: process with  1  covariance structure(s)
#> grf: nugget effect is: tausq= 0 
#> grf: covariance model 1 is a pure nugget effect
#> grf: decomposition algorithm used is:  cholesky 
#> grf: End of simulation procedure. Number of realizations: 1
points.geodata(sim3, main = "Independencia", col = tim.colors(), cex.max = 3)

Comentar

par(mfrow = c(1, 2))
set.seed(2022)
sim2 <- grf(441, grid = "reg", cov.pars = c(1, 0.25))
#> grf: generating grid  21  *  21  with  441  points
#> grf: process with  1  covariance structure(s)
#> grf: nugget effect is: tausq= 0 
#> grf: covariance model 1 is: exponential(sigmasq=1, phi=0.25)
#> grf: decomposition algorithm used is:  cholesky 
#> grf: End of simulation procedure. Number of realizations: 1
image(sim2, main = "Dependencia espacial (positiva) (ii)", col = tim.colors())


# Independencia
set.seed(2022)
sim4 <- grf(441, grid = "reg", cov.pars = c(0.01, 0))
#> grf: generating grid  21  *  21  with  441  points
#> grf: process with  1  covariance structure(s)
#> grf: nugget effect is: tausq= 0 
#> grf: covariance model 1 is a pure nugget effect
#> grf: decomposition algorithm used is:  cholesky 
#> grf: End of simulation procedure. Number of realizations: 1
image(sim4, main = "Independencia", col = tim.colors())

3.2 Datos espaciales GEMA

Los datos espaciales, también conocidos como datos geoespaciales, son aquellos datos relacionados o que contienen información de una localización o área geográfica de la superficie de la Tierra.

La forma más intuitiva de representar los datos espaciales es a través de un mapa.

<!-- Propuesta: mapa temático cualquiera con poco código para empezar como el libro de SDS  https://keen-swartz-3146c4.netlify.app/intro.html#a-first-map  -->
<!--    Y comentar algunas características de los datos espaciales que luego se explican. -->

# Mapa de porcentaje de mujeres en Castilla-La Mancha

library(mapSpain)

# Datos de población
pob <- mapSpain::pobmun19

# Porcentaje
pob$porc_mujeres <- pob$women / pob$pob19

# Datos geoespaciales
geo <- esp_get_munic(region = "Castilla-La Mancha")

# Une ambos datos
geo_pob <- merge(geo,
  pob,
  by = c("cpro", "cmun"),
  all.x = TRUE
)

# Mapa básico
plot(geo_pob["porc_mujeres"],
  # Cambiamos titulo
  main = "Castilla-La Mancha: % mujeres (2019)",

  # Cambiamos la paleta de colores para hacerlo mas atractivo
  border = NA,
  pal = hcl.colors(12, "RdYlBu")
)

3.3 Clasificación de datos espaciales

Tal y como acabamos de señalar y de acuerdo con Schabenberger y Gotway (2005, p. 6), debido a que los datos espaciales surgen en una gran variedad de campos y aplicaciones, también hay una gran variedad de tipos de datos espaciales, estructuras y escenarios. Por tanto, una clasificación exhaustiva de los datos espaciales sería un reto muy difícil y hemos apostado por una clasificación general, simple y útil de datos espaciales proporcionada por Cressie (1993).

La clasificación de Cressie de datos espaciales se basa en la naturaleza del dominio espacial en estudio. Dependiendo de esto, podemos tener: datos geoestadísticos, datos de patrones de puntos y datos latice.

Siguiendo a Cressie (1993), sea $s ∈ ℝ^d$ una localización en un espacio Euclideo \(d-\)dimensional y ${Z(s)∶ s ∈ ℝ^d}$ una función aleatoria espacial, donde \(Z\) representa el atributo en el cual estamos interesados:

1. Datos geoestadísticos: Surgen cuando el dominio en estudio es conjunto y fijo \(D\). Es decir: (i) \(Z(s)\) se puede observar en cualquier punto del dominio (continuo); y (ii) los puntos en \(D\) no son estocásticos (son fijos, \(D\) es el mismo para todas las realizaciones de la función aleatoria espacial ).

   Algunos ejemplos de datos geoestadísticos son el nivel de un contaminante en una ciudad, los valores de precipitación o temperatura del aire en un país, las concentraciones de metales pesados en la capa superior del suelo de una región, etc.

   Es obvio que, al menos en teoría, el nivel de un contaminante específico podría medirse en cualquier lugar de la ciudad; Lo mismo puede decirse de las mediciones de precipitaciones o temperaturas del aire en un país o concentraciones de un metal pesado en una región.

   Sin embargo, en la práctica, no es posible una observación exhaustiva del proceso espacial. Por lo general, el proceso espacial se observa en un conjunto de ubicaciones (por ejemplo, el nivel de un contaminante específico en una ciudad se observa en los puntos donde están ubicadas las estaciones de monitoreo) y, basado en tales valores observados, el análisis geoestadístico reproduce el comportamiento de el proceso espacial en todo el dominio de interés.

   En el análisis geoestadístico lo más importante es cuantificar la correlación espacial entre observaciones (a través de la herramienta básica en geoestadística, el semivariograma) y utilizar esta información para lograr los objetivos anteriores.

# ejemplo

2. Datos reticulares: Surgen cuando: (i) el dominio bajo estudio \(D\) es discreto, es decir, \(Z(s)\) puede observarse en una serie de ubicaciones fijas que pueden enumerarse. Estas ubicaciones pueden ser puntos o regiones, pero generalmente son códigos postales, pistas censales, vecindarios, provincias, países, etc., y los datos en la mayoría de los casos son datos agregados espacialmente sobre estas áreas. Aunque estas regiones pueden tener una forma regular, normalmente la forma que tienen es irregular, y esto, junto con el carácter espacialmente agregado de la datos, es por lo que los datos latice tambien se denominan datos regionales. Y (ii) las ubicaciones en \(D\) no son estocásticas. Por supuesto, un concepto clave en el análisis de los datos lattice es el vecindario.

   Algunos ejemplos de reticulares incluyen la tasa de desempleo por estados, los datos de delincuencia por comarcas, rendimientos agrícolas en parcelas, precios medios de la vivienda por provincias, etc.

# ejemplo

3. Procesos de puntos: Mientras que en los datos geoestadísticos y reticulares el dominio \(D\) es fijo, en los datos de patrones puntuales el dominio es discreto o continuo, pero aleatorio. Los patrones de puntos surgen cuando el atributo bajo estudio es la ubicación de los eventos (observaciones). Es decir, el interés radica en dónde ocurren eventos de interés.

   Algunos ejemplos de patrones de puntos son la ubicación de incendios en una región española, la ubicación de los árboles en un bosque o la ubicación de nidos en una colonia de aves reproductoras, la localización de los delitos en una ciudad, entre muchas otras.

   En estos En los casos, es obvio que D es aleatorio y los puntos de observación no dependen del investigador. El principal objetivo del análisis de patrones de puntos es determinar si la ubicación de los eventos tiende a exhibir un patrón sistemático sobre el área en estudio o, por el contrario, son aleatoriamente repartido.

   Más concretamente, nos interesa analizar si la ubicación de los eventos es completamente aleatorio espacialmente (la ubicación donde ocurren los eventos no se ve afectada por la ubicación de otros eventos), uniforme o regular (cada punto está tan lejos de todos sus vecinos como sea posible) o agrupados o agregados (la ubicación de los eventos se concentra en grupos).

# ejemplo

Los formatos más comunes de datos espaciales son vectores y ráster.

Diego

4 Formatos de datos espaciales IGUAL SE PUEDE PONER UN NOMBRE MÁS GENÉRICO Y QUE VAYA TODOD

4.1 Vectores

4.2 Ráster.

4.3 CRS

Un sistema de referencia de coordenadas (o CRS por sus siglas en inglés, Coordinate Reference System) permite relacionar datos espaciales con su localización en la superficie terrestre. Los CRS constituyen por tanto un aspecto fundamental en el análisis y representación de datos espaciales, ya que nos permiten identificar con exactitud la posición de los datos sobre el globo terráqueo.

Así mismo, cuando se trabaja con datos espaciales provenientes de distintas fuentes de información, es necesario comprobar que dichos datos se encuentran definidos en el mismo CRS:

En el ejemplo anterior, ambos puntos (verde y rojo) presentan los mismos valores de coordenadas en los ejes X e Y, en este caso las correspondientes a la ciudad de Toledo.

Sin embargo, presentan distintos CRS. Por este motivo, al representar ambos puntos en un mapa, se observa que no se están refiriendo a la misma localización geográfica. Esto es así porque el CRS define la referencia (punto x=0 e y =0) y las unidades de los ejes (grados, metros, millas).

Como conclusión, además de disponer de las coordenadas de los datos espaciales, es necesario conocer el CRS en el que están definidos para conocer de manera exacta su localización geográfica. Además, nótese que para cualquier análisis de datos espaciales es necesario que todos los geodatos se encuentren referenciados en el mismo CRS. Esto se consigue transformando (o proyectando) los datos a un CRS común, nunca sobreescribiendo el CRS de los mismos.

4.3.1 Tipos de CRS

A continuación se definen los dos grandes tipos de CRS, los CRS geográficos y los CRS proyectados.

4.3.1.1 CRS geográficos

Los CRS geográficos son aquellos en los que los parámetros empleados para localizar una posición espacial son la latitud y la longitud:

• Latitud: Es la distancia angular expresada en grados sobre el plano definido por el ecuador terrestre. Determina la posición sobre de una localización en el eje Norte-Sur de la Tierra y toma valores en el rango \([-90º,90º]\) . Las líneas imaginarias determinadas por una sucesión de puntos con la misma latitud a lo largo del eje Este-Oeste se denominan paralelos:

• Longitud: Es la distancia angular expresada en grados sobre el plano definido por el meridiano de Greenwich. Determina la posición sobre de una localización en el eje Este-Oeste de la Tierra y toma valores en el rango \([-180º,180º]\) . Las líneas imaginarias determinadas por una sucesión de puntos con la misma longitud a lo largo del eje Este-Oeste se denominan meridianos:

Es muy importante destacar que en un sistema de coordenadas geográfico, es decir, basado en latitudes y longitudes, las distancias entre dos puntos representan distancias angulares. Por ejemplo, la distancia entre el meridiano de Greenwich y el meridiano correspondiente a la longitud 20º siempre es de +20º. Sin embargo, debido a la forma esférica de la Tierra, la longitud en metros entre ambos meridianos no es constante:

4.3.1.2 CRS proyectados

La representación de formas tridimensionales en un soporte plano (dos dimensiones) presenta algunos retos. Por ello, es habitual trabajar con proyecciones de mapas.

Una proyección geográfica es un método para reducir la superficie de la esfera terrestre a un sistema cartesiano de dos dimensiones. Para ello, es necesario transformar las coordenadas longitud y latitud en coordenadas cartesianas x e y.

Es importante destacar que las proyecciones pueden incluir un punto de origen (X=0, Y=0) y unas unidades de distancia (habitualmente metros) específicas. Por ejemplo, la proyección cónica equiáreas de Albers (específica para Estados Unidos) define su punto de referencia (0,0) en la latitud 40º N y longitud 96º, y la unidad de variación están definida en metros. De ahí la importancia de conocer el CRS de los datos geográficos, como se expuso al principio de este tema.

Existen varias familias de proyecciones, que se pueden clasificar de diversas maneras:

4.3.1.2.1 Por tipo de superficie de proyección

El proceso de trasladar puntos de una esfera a un plano puede plantearse de manera práctica como el ejercicio de envolver una esfera con una superificie plana (como una hoja de papel) y trasladar los puntos de la esfera de manera lineal al punto de la superficie plana más cercano a ella.

A partir de este ejercicio, se plantean tres posibles soluciones, dependiendo del tipo de superficie que se use para proyectar:

• Proyecciones cilíndricas: Son aquellas proyecciones donde la superficie de proyección conforma un cilindro alrededor de la Tierra. Una de las proyecciones cilíndricas más conocidas es la proyección de Mercator.

• Proyecciones cónicas: En este tipo de proyecciones, se plantea la superficie de proyección como una forma cónica. Como ejemplo, la proyección cónica equiáreas de Albers es una de las proyecciones que más suele usarse en la representación de mapas de América del Norte:

  

• Proyecciones acimutales o planares: En este tipo de proyección se proyecta una porción de la Tierra sobre un plano que es tangente a la misma en el punto de referencia. Como ejemplos de proyecciones acimutales podemos destacar la proyección ortográfica:

  

4.3.1.2.2 Por métrica a preservar

Es importante tener en cuenta que cualquier proyección de la superficie de la Tierra produce distorsiones en una o varias características geográficas. Como ejemplos clásicos, la proyección de Mercator produce distorsiones del área especialmente en aquellas regiones más cercanas a los polos (Groenlandia, que la proyección de Mercator presenta una área similar a la de África, presenta menor superificie real que Argelia). Otras de las métricas que suele verse distorsionada son la distancia entre dos puntos geográficos, la dirección o la forma de regiones de la Tierra.

A lo largo de la Historia se han desarrollado diversas proyecciones cuyo objetivo es preservar alguna o varias de las propiedades mencionadas anteriormente, sin embargo es importante destacar que no existe una proyección que sea capaz de preservar todas las métricas a la vez.

Según la metrica a presevar, las proyecciones se pueden clasificar en:

• Proyecciones conformales: Estas proyecciones intentan preservar los ángulos que se forman en la superficie terrestre. Por ejemplo, la proyección de Mercator representa ángulos rectos en las intersecciones de los paralelos y los meridianos.

  

• Proyecciones equivalentes: Estas proyecciones preservan las proporciones de las áreas, provocando a su vez deformaciones en el resto de características, como la forma o los ángulos. La proyección acimutal equivalente de Lambers es un tipo de proyección equivalente.

  

• Proyecciones equidistantes: Este tipo de proyección preserva la distancia entre dos puntos geográficos específicos. Por ejemplo, la proyección Plate carré preserva la distancia entre el Polo Norte y el Polo Sur.

• Proyecciones de compromiso: Este tipo de proyección no intenta preservar ninguna métrica en concreto. En su lugar, se centran en intentar encontrar un equilibrio entre las diversas distorsiones que provocan para intentar dar una representación más o menos representativa de la superficie terrestre. La proyección de Winkel Tripel, usada en los mapas de National Geographic, es un ejemplo de proyección de compromiso.

En los anteriores ejemplos se ha añadido a cada proyección la indicatriz de Tissot. Consiste en una serie de círculos imaginarios de igual área distribuidos sobre la superficie esférica de la Tierra en determinados puntos. De este manera, al presentar la indicatriz de Tissot en una proyección específica, se puede entender de una manera intuitiva la distorsión provocada por dicha proyección, ya que los círculos se ven distorsionados o preservados según los parámetros y la naturaleza de la proyección en cuestión.

4.3.2 Trabajando con proyecciones en R

Existe toda una serie de proyecciones predefinidas, identificadas mediante los códigos EPSG, ESRI, WKT o proj4 (en desuso en R, pero todavía admitidos). Existen varios recursos web donde se pueden consultar y seleccionar los códigos correspondientes:

• https://epsg.io/

• https://spatialreference.org/

• https://proj.org/operations/projections/index.html

El paquete sf permite obtener los parámetros de estas proyecciones mediante la función st_crs():

library(sf)

# Ejemplo: EPSG WGS 84 (Sistema Global GPS): EPSG 4326

st_crs(4326)
#> Coordinate Reference System:
#>   User input: EPSG:4326 
#>   wkt:
#> GEOGCRS["WGS 84",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["geodetic latitude (Lat)",north,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         AXIS["geodetic longitude (Lon)",east,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>     USAGE[
#>         SCOPE["unknown"],
#>         AREA["World"],
#>         BBOX[-90,-180,90,180]],
#>     ID["EPSG",4326]]

# Usando código ESRI North America Albers Equal Area Conic

st_crs("ESRI:102008")
#> Coordinate Reference System:
#>   User input: ESRI:102008 
#>   wkt:
#> PROJCRS["North_America_Albers_Equal_Area_Conic",
#>     BASEGEOGCRS["NAD83",
#>         DATUM["North American Datum 1983",
#>             ELLIPSOID["GRS 1980",6378137,298.257222101,
#>                 LENGTHUNIT["metre",1]]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["Degree",0.0174532925199433]]],
#>     CONVERSION["North_America_Albers_Equal_Area_Conic",
#>         METHOD["Albers Equal Area",
#>             ID["EPSG",9822]],
#>         PARAMETER["Latitude of false origin",40,
#>             ANGLEUNIT["Degree",0.0174532925199433],
#>             ID["EPSG",8821]],
#>         PARAMETER["Longitude of false origin",-96,
#>             ANGLEUNIT["Degree",0.0174532925199433],
#>             ID["EPSG",8822]],
#>         PARAMETER["Latitude of 1st standard parallel",20,
#>             ANGLEUNIT["Degree",0.0174532925199433],
#>             ID["EPSG",8823]],
#>         PARAMETER["Latitude of 2nd standard parallel",60,
#>             ANGLEUNIT["Degree",0.0174532925199433],
#>             ID["EPSG",8824]],
#>         PARAMETER["Easting at false origin",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8826]],
#>         PARAMETER["Northing at false origin",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8827]]],
#>     CS[Cartesian,2],
#>         AXIS["(E)",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1]],
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1]],
#>     USAGE[
#>         SCOPE["unknown"],
#>         AREA["North America - Canada and USA (CONUS, Alaska mainland)"],
#>         BBOX[23.81,-172.54,86.46,-47.74]],
#>     ID["ESRI",102008]]

# Usando proj4string: Robinson

st_crs("+proj=robin")
#> Coordinate Reference System:
#>   User input: +proj=robin 
#>   wkt:
#> PROJCRS["unknown",
#>     BASEGEOGCRS["unknown",
#>         DATUM["World Geodetic System 1984",
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1]],
#>             ID["EPSG",6326]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8901]]],
#>     CONVERSION["unknown",
#>         METHOD["Robinson"],
#>         PARAMETER["Longitude of natural origin",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8802]],
#>         PARAMETER["False easting",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8806]],
#>         PARAMETER["False northing",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8807]]],
#>     CS[Cartesian,2],
#>         AXIS["(E)",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]],
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]]]

De esta manera, es posible proyectar un objeto sf mediante la función st_transform():

# Usa datos del paquete mapSpain

library(giscoR)

paises <- gisco_get_countries()

# Comprobamos el CRS de estos datos
# Se puede almacenar en un objeto y usar posteriormente
st_crs(paises)
#> Coordinate Reference System:
#>   User input: EPSG:4326 
#>   wkt:
#> GEOGCS["WGS 84",
#>     DATUM["WGS_1984",
#>         SPHEROID["WGS 84",6378137,298.257223563,
#>             AUTHORITY["EPSG","7030"]],
#>         AUTHORITY["EPSG","6326"]],
#>     PRIMEM["Greenwich",0,
#>         AUTHORITY["EPSG","8901"]],
#>     UNIT["degree",0.0174532925199433,
#>         AUTHORITY["EPSG","9122"]],
#>     AUTHORITY["EPSG","4326"]]

# Plot
plot(st_geometry(paises), axes = TRUE)

# Proyectamos a Mercator
# El eje cambia porque Mercator usa metros
paises_merc <- st_transform(paises, st_crs(3857))
plot(st_geometry(paises_merc), axes = TRUE)

# Proyectamos a Robinson
paises_robin <- st_transform(paises, st_crs("+proj=robin"))
plot(st_geometry(paises_robin), axes = TRUE)

Como se comentó anteriormente, cuando se usan geodatos de diversas fuentes, es necesario que todos presenten el mismo CRS. En este ejemplo se muestra lo que ocurre si esto no se cumple:

# Añadimos a este mapa puertos mundiales de giscoR

puertos <- gisco_get_ports()
plot(st_geometry(paises_robin), main = "Puertos en el mundo")
plot(st_geometry(puertos), add = TRUE, col = "red", pch = 20)

# Ha habido algun error... Comprueba CRS

st_crs(puertos) == st_crs(paises_robin)
#> [1] FALSE

# Los puertos no están en Robinson! Proyectamos al mismo CRS
puertos_robin <- st_transform(puertos, st_crs(paises_robin))
plot(st_geometry(paises_robin), main = "Puertos en el mundo")
plot(st_geometry(puertos_robin), add = TRUE, col = "blue", pch = 20)

Como vemos, en el primer mapa los puertos se concentran en un único punto, dado que no están referenciados en el mismo CRS. Tras proyectarlos, el mapa se representa adecuadamente.

En otros paquetes, como sp o raster, existen funciones parecidas. Cuando empleemos el paquete sp podemos usar las funciones CRS() y spTransform():

library(sp)

# Convertimos sf a sp
paises_sp <- as(paises, "Spatial")
paises_sp@proj4string
#> Coordinate Reference System:
#> Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
#> WKT2 2019 representation:
#> GEOGCS["WGS 84",
#>     DATUM["WGS_1984",
#>         SPHEROID["WGS 84",6378137,298.257223563,
#>             AUTHORITY["EPSG","7030"]],
#>         AUTHORITY["EPSG","6326"]],
#>     PRIMEM["Greenwich",0,
#>         AUTHORITY["EPSG","8901"]],
#>     UNIT["degree",0.0174532925199433,
#>         AUTHORITY["EPSG","9122"]],
#>     AUTHORITY["EPSG","4326"]]

# En sp podemos usar:
CRS("+proj=robin")
#> Coordinate Reference System:
#> Deprecated Proj.4 representation:
#>  +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#> WKT2 2019 representation:
#> PROJCRS["unknown",
#>     BASEGEOGCRS["unknown",
#>         DATUM["World Geodetic System 1984",
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1]],
#>             ID["EPSG",6326]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8901]]],
#>     CONVERSION["unknown",
#>         METHOD["Robinson"],
#>         PARAMETER["Longitude of natural origin",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8802]],
#>         PARAMETER["False easting",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8806]],
#>         PARAMETER["False northing",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8807]]],
#>     CS[Cartesian,2],
#>         AXIS["(E)",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]],
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]]]

# O también desde sf
CRS(st_crs(paises_robin)$proj4string)
#> Coordinate Reference System:
#> Deprecated Proj.4 representation:
#>  +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#> WKT2 2019 representation:
#> PROJCRS["unknown",
#>     BASEGEOGCRS["unknown",
#>         DATUM["World Geodetic System 1984",
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1]],
#>             ID["EPSG",6326]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8901]]],
#>     CONVERSION["unknown",
#>         METHOD["Robinson"],
#>         PARAMETER["Longitude of natural origin",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8802]],
#>         PARAMETER["False easting",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8806]],
#>         PARAMETER["False northing",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8807]]],
#>     CS[Cartesian,2],
#>         AXIS["(E)",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]],
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]]]


paises_sp_robin <- spTransform(paises_sp, CRS("+proj=robin"))
plot(paises_sp_robin)

En el caso de un objeto raster, podemos usar crs() y projectRaster():

library(raster)


# Extrae información de altitud para Paises Bajos
elev <- getData("alt", country = "NLD", path = tempdir())

raster::crs(elev)
#> Coordinate Reference System:
#> Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
#> WKT2 2019 representation:
#> GEOGCRS["unknown",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]],
#>         ID["EPSG",6326]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433],
#>         ID["EPSG",8901]],
#>     CS[ellipsoidal,2],
#>         AXIS["longitude",east,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]],
#>         AXIS["latitude",north,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]]]

plot(elev)

# Transforma

elev_robinson <- projectRaster(elev, crs = crs("+proj=robin"))
plot(elev_robinson)

Por último, en el paquete terra las funciones correspondientes son crs() y project():

library(terra)

# Convierte de raster a terra
elev_terra <- rast(elev)
cat(terra::crs(elev_terra))
#> GEOGCRS["unknown",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]],
#>         ID["EPSG",6326]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433],
#>         ID["EPSG",8901]],
#>     CS[ellipsoidal,2],
#>         AXIS["longitude",east,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]],
#>         AXIS["latitude",north,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]]]

plot(elev_terra)

# Transforma
elev_terra_robinson <- terra::project(elev_terra, terra::crs(elev_terra))
plot(elev_terra_robinson)

4.3.3 ¿Qué proyección uso?

El CRS adecuado para cada análisis depende de la localización y el rango espacial de los datos. Un CRS adecuado para representar un mapa del mundo puede no serlo para representar datos de zonas específicas de la Tierra. Los recursos web mencionados anteriormente permiten la búsqueda de CRS por zona geográfica, y adicionalmente en R existe el paquete crsuggest (Walker, 2021) que nos facilita la labor, sugiriendo el CRS más adecuado para cada zona:

library(crsuggest)

# CRS para Paises Bajos

# Usando raster
suggest_crs(elev)
#> # A tibble: 10 × 6
#>    crs_code crs_name         crs_type crs_gcs crs_units crs_proj4               
#>    <chr>    <chr>            <chr>      <dbl> <chr>     <chr>                   
#>  1 28992    Amersfoort / RD… project…    4289 m         +proj=sterea +lat_0=52.…
#>  2 28991    Amersfoort / RD… project…    4289 m         +proj=sterea +lat_0=52.…
#>  3 5651     ETRS89 / UTM zo… project…    4258 m         +proj=tmerc +lat_0=0 +l…
#>  4 5649     ETRS89 / UTM zo… project…    4258 m         +proj=tmerc +lat_0=0 +l…
#>  5 3812     ETRS89 / Belgia… project…    4258 m         +proj=lcc +lat_0=50.797…
#>  6 3447     ETRS89 / Belgia… project…    4258 m         +proj=lcc +lat_0=50.797…
#>  7 31370    Belge 1972 / Be… project…    4313 m         +proj=lcc +lat_0=90 +lo…
#>  8 31300    Belge 1972 / Be… project…    4313 m         +proj=lcc +lat_0=90 +lo…
#>  9 21500    Belge 1950 (Bru… project…    4809 m         +proj=lcc +lat_0=90 +lo…
#> 10 23095    ED50 / TM 5 NE   project…    4230 m         +proj=tmerc +lat_0=0 +l…

# Probamos sugerencia
crs_suggest <- suggest_crs(elev, limit = 1)

elev_suggest <- projectRaster(elev, crs = raster::crs(crs_suggest$crs_proj4))

plot(elev_suggest)

# Ejemplo con sf: China

china <- gisco_get_countries(country = "China")
plot(st_geometry(china), axes = TRUE)

china_crs <- suggest_crs(china, limit = 1)


china_suggest <- st_transform(
  china,
  st_crs(as.integer(china_crs$crs_code))
)


plot(st_geometry(china_suggest), axes = TRUE)

4.4 Tipos de ficheros (ver contigo).

5 APLICACIÓNES:

para explicar cómo se hace un merge y luego se representa en un mapa o en unos puntos

6 EXTENSIONES:

ESPACIO-TEMPORAL SÍ. Pensar qué. Temperatura mínima? Mapas de contorno y en 3D, se usan mucho.

Leaflet: el de contaminación con las capas quedó muy bien, se puede poner.

bbox, elecciones (discreta), google satélite ? SE podía poner algo de COVID…

Todo muy sencillo y con un ejemplito

6.1 Mapas temáticos (esto podría ir dentro de los ejemplos…, ver con DIEGO)

mapSpain, tmap Datos que ya tenemos: -Medioambiente Madrid. -Temperatura mínima ESpaña - Alguno económico a nivel municipal para españa (PIB, PARO,..) - Alguno de La Palma.

References

Abrahart, R. J., Openshaw, S., Abrahart, R. J., & See, L. M. (Eds.). (2000). Geocomputation. CRC Press. https://doi.org/10.4324/9780203305805

Cressie, N. A. C. (1993). Statistics for spatial data. John Wiley & Sons, Inc. https://doi.org/10.1002/9781119115151

Lovelace, R., Nowosad, J., & Muenchow, J. (2019). Geocomputation with R. CRC Press.

Walker, K. (2021). crsuggest: Obtain suggested coordinate reference system information for spatial data.